Matlab Codes

Chapter 4

Matlab code 4.3: Matlab file “Figure 4-8.m”

%--------------------------------------------------------------------

% This code can be used to generate Figure 4.8

%--------------------------------------------------------------------

 

clear all;

close all;

clc

%% The comparison of m-D curve obtained by monostatic radar and bistatic radar

c = 3e8;

j = sqrt(-1);

fc = 10e9; % carrier frequency of transmitted signal

radt = [0 0 0]; % coordinates of transmitter

radr = [7000 0 0]; % coordinates of receiver

cord = 1000*[3 4 5]; % coordinates of target centre in the world coordinate system

colo = [0 0 0]; % local coordinate system's origin in the world coordinate system

R0 = colo-cord; % distance between target and radar(monostatic radar)

Rt = radt-cord; % distance between target and transmitter (bistatic radar)

Rr = cord-radr; % distance between target and receiver (bistatic radar)

w = pi*[1 2 1]; % angular velocity of rotation

tp = 4e-6; % pulse duration

B = 500e6; % bandwidth

r = [2 2 -6]; % rotation radius vector of the scaterers

t = 1.5; % radar illumimated time

v = 0; % translational velocity of target

prf = 1000; % pulse repetition frequency

pri = 1/prf; % pulse repetition interval

dt = 0:pri:t-pri; % slow time sampling interval

m = length(dt); % number of time sampling

ae = [0 pi/4 pi/5]; % Euler angle

ri1 = [cos(ae(3)) sin(ae(3)) 0;-sin(ae(3)) cos(ae(3)) 0;0 0 1];

ri2 = [1 0 0;0 cos(ae(2)) sin(ae(2));0 -sin(ae(2)) cos(ae(2))];

ri3 = [cos(ae(1)) sin(ae(1)) 0;-sin(ae(1)) cos(ae(1)) 0;0 0 1];

ri = ri1*ri2*ri3; % initial rotation matrix

r0 = r-colo; % scatterers in the local coordinate system

r0r = ri*r0'; % scatterers in the reference coordinate system

w = ri*w'; % angular velocity of rotation in the reference coordinate system

omega = sqrt(sum(w.^2)); % angular frequency

we = w/omega; % unit vector of rotation

wr = [0 -we(3) we(2);we(3) 0 -we(1);-we(2) we(1) 0]; % skew symmetric matrix

Rr1 = zeros(length(r),length(dt)); % distance between the scatterers and radar (monostatic radar)

Rr2 = zeros(length(r),length(dt)); % distance between the scatterers and radar (bistatic radar)

for i = 1:m

    % the instantaneous range in monostatic radar

    Rr1(:,i) = sqrt((R0'+v*dt(i)+((eye(3)+wr*sin(omega*dt(i))+wr^2*(1-cos(omega*dt(i))))*r0r))'...

                   *(R0'+v*dt(i)+((eye(3)+wr*sin(omega*dt(i))+wr^2*(1-cos(omega*dt(i))))*r0r)));

    % the instantaneous range in bistatic radar

    Rr2(:,i) = sqrt((Rt'+v*dt(i)+((eye(3)+wr*sin(omega*dt(i))+wr^2*(1-cos(omega*dt(i))))*r0r))'...

                   *(Rt'+v*dt(i)+((eye(3)+wr*sin(omega*dt(i))+wr^2*(1-cos(omega*dt(i))))*r0r)))...

              +sqrt((Rr'+v*dt(i)+((eye(3)+wr*sin(omega*dt(i))+wr^2*(1-cos(omega*dt(i))))*r0r))'...

                   *(Rr'+v*dt(i)+((eye(3)+wr*sin(omega*dt(i))+wr^2*(1-cos(omega*dt(i))))*r0r)));

end

%% monostatic radar

rmin = sqrt(sum((cord-colo).^2)); % the minmum detectable distance

rmax = sqrt(sum((cord-colo).^2))+sqrt(sum((r-colo).^2)); % the maximum detectable distance

ts = rmin/c; % the starting point of the time sampling

te = rmax/c+tp; % the ending point of the time sampling

mu = B/tp; % frequency modulation rate

fs = B; % fast time sampling frequency

kdt = 1/fs; % fast time sampling interval

nt = ceil((te-ts)/(2*kdt)); % number of fast time sampling

df = fs/nt; % number of frequency sampling

tk = ts+(0:nt-1)*kdt; % fast time

s1 = zeros(nt,length(dt)); % echo signal matrix (monostatic radar)

for n = 1:length(r)

    td1 = tk(:)*ones(1,m)-2*ones(nt,1)*Rr1(n,:)/c; % echo delay

    s1 = s1+exp(j*2*pi*fc*td1+j*pi*mu*td1.^2);

end

Rr0 = sqrt(sum(R0.^2)); % reference range

td0 = tk(:)*ones(1,m)-2*Rr0/c; % time delay of reference range

sref = exp(j*2*pi*fc*td0+j*pi*mu*td0.^2); % reference signal

sdt1 = s1.*conj(sref); % dechirp

sdf1 = fftshift(fft(sdt1),1);   

clear td1 tk

 

x = dt;

y = -(df*(0:nt-1)-B/2)*c/(2*mu); % turn the fast time frequency-slow time plane to the range-slow time plane

figure(1)

contour(x,y,abs(sdf1),5)

xlabel('Slow time (s)')

ylabel('Range (m)')

axis([0,1.5,-10,10])

clear td0 sref s1 sdt1 sdf1 x y

 

%% bistatic radar

Width = 10; % the imaging scene width

rmin_1 = sqrt(sum((cord-colo).^2))+sqrt(sum((cord-radr).^2))-Width; % the minmum detectable distance

rmax_1 = sqrt(sum((cord-colo).^2))+sqrt(sum((cord-radr).^2))+Width; % the maximum detectable distance

ts_1 = rmin_1/c; % the starting point of the time sampling

te_1 = rmax_1/c+tp; % the ending point of the time sampling

mu = B/tp; % frequency modulation rate

fs = B; % fast time sampling frequency

kdt = 1/fs; % fast time sampling interval

nt = 2*ceil((te_1-ts_1)/(2*kdt)); % number of fast time sampling

df = fs/nt; % number of frequency sampling

tk_1 = ts_1+(0:nt-1)*kdt; % fast time

s2 = zeros(nt,length(dt)); % echo signal matrix (bistatic radar)

for n = 1:length(r)

    td2 = tk_1(:)*ones(1,m)-ones(nt,1)*Rr2(n,:)/c; % echo delay

    s2 = s2+exp(j*2*pi*fc*td2+j*pi*mu*td2.^2);

end

clear td2

Rr1 = sqrt(sum(R0.^2))+sqrt(sum(Rr.^2)); % reference range

td1 = tk_1(:)*ones(1,m)-Rr1/c;

sref1 = exp(j*2*pi*fc*td1+j*pi*mu*td1.^2); % reference signal

sdt2 = s2.*conj(sref1); % dechirp

sdf2 = fftshift(fft(sdt2),1);

clear td1 sref1 sdt2

x = dt;

y = -(df*(0:nt-1)-B/2)*c/mu;

figure(2)

contour(x,y,abs(sdf2),5)

xlabel('Slow time (s)')

ylabel('Range (m)')

axis([0,1.5,-10,10])